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ABSTRACT 


The  track  keeping  characteristics  of  an  autonomous  ocean  vehicle  in  the  presence 
of  a  realistic  time  lag  on  the  vehicle’s  positional  information  is  examined  in  two  ways. 
A  Hopf  bifurcation  analysis  is  applied  to  better  predict  vehicle  behavior  within  its  region 
of  linearized  classical  stability.  Additionally,  the  effects  on  guidance/ control  stability  by 
the  use  of  a  single  stage  memory  model  incorporating  the  two  previous  position  data 
points  in  place  of  a  single  position  data  point  is  investigated.  Results  are  presented  using 
a  dynamic  model  of  the  Naval  Postgraduate  School  Autonomous  Underwater  Vehicle  II 
(NPS  AUV  D). 
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I. 


INTRODUCTION 


A.  BACKGROUND 

Autonomous  vehicles  suitable  for  use  in  modern 
applications  require  high  maneuverability ,  quick  response,  and 
robust  performance  characteristics  [Ref.  7].  In  order  to 
maintain  accurate  track  keeping  characteristics,  a  suitable 
combination  of  path  planning,  navigation,  guidance,  and 
control  design  is  needed  [Ref.  10]  .  Sufficient  information  is 
obtained  from  charted  obstacles  and  the  environment  for  smooth 
paths  to  be  generated  for  the  vehicle  to  follow  [Ref.  8]. 
Although  it  is  possible  to  design  a  nonlinear  controller  which 
incorporates  and  utilizes  information  on  the  nonlinear  dynamic 
properties  of  the  vehicle,  as  well  as  the  geometric 
nonlinearities  of  the  desired  track,  the  resulting  scheme 
tends  to  be  very  complex  and  time  consuming  [Ref.  3].  The 
alternative  is  to  use  separated  navigation,  guidance,  and 
autopilot  functions.  A  certain  level  of  feedback  is  provided 
at  the  navigation  level  through  the  use  of  sonars  in  order  to 
replan  a  path  due  to  uncharted  obstacles  or  changes  in  mission 
objectives.  This  operation  is  event -driven  and  occurs 
asynchronously,  and  in  many  cases  it  does  not  affect  stability 
and  performance  of  the  overall  navigation,  guidance,  and 
control  scheme.  Based  on  the  provided  navigational 


information,  the  guidance  law  generates  heading  commands, 
which  are  executed  by  the  control  law  by  appropriate  use  of 
vehicle  effectors,  such  as  control  surfaces  and  cross  body 
thrusters.  However,  the  time  required  to  process  sonar  data 
and/or  inertial  position  information  may  be  significant  and 
cannot  be  neglected  [Ref.  9].  In  addition,  the  guidance  and 
control  laws  must  be  as  fast  as  possible  in  order  to  ensure 
accurate  path  keeping  characteristics.  The  navigational 
position  information  time  lags,  as  well  as  the  dynamic  lags 
due  to  the  vehicle  inertia,  however,  set  a  limit  on  the 
vehicle  reaction  time.  Therefore,  stability  of  the  combined 
scheme  becomes  an  issue  which  requires  analysis. 

B.  OBJECTIVE 

Previous  studies  have  established  boundaries  for  guidance 
and  control  laws  in  the  horizontal  [Ref.  10]  and  vertical 
plane  [Ref.  12]  maneuvering  along  straight  line  paths,  as  well 
as  curved  segments  [Ref.  11].  The  most  important  assumption 
in  those  results  was  the  existence  of  instantaneous  positional 
information  updates  when  needed.  In  this  work  we  relax  the 
latter  assumption.  Stability  boundaries  are  computed 
parametrized  by  the  amount  of  positional  information  time  lag. 
Results  are  presented  based  on  eigenvalue  and  frequency 
response  techniques.  This  thesis  builds  on  the  linear 
stability  analysis  performed  in  [Ref.  13]  which  can  predict 
the  shape  of  the  stability  boundaries.  In  this  work,  vehicle 


motions  in  the  vicinity  of  the  corresponding  boundary  are 
assessed  using  nonlinear  bifurcation  theory  [Ref.  5].  It  is 
shown  that  the  loss  of  stability  occurs  as  generic  Hopf 
bifurcations,  where  upon  loss  of  stability  of  straight  line 
motion  a  family  of  periodic  solutions  appears.  Nonlinear 
expansions  utilizing  center  manifold  approximations  and 
integral  averaging  techniques,  reveal  that  these  Hopf 
bifurcations  can  occur  either  in  supercritical  or  subcritical 
forms.  In  the  supercritical  case,  a  stable  family  of  limit 
cycles  is  generated  immediately  after  the  nominal  motion 
becomes  unstable.  In  the  subcritical  case,  however,  the 
resulting  limit  cycles  are  initially  unstable  and  they  are 
generated  even  before  the  nominal  motion  loses  its  stability. 
In  the  latter  case,  the  domain  of  convergence  of  straight  line 
motion  becomes  increasingly  smaller  as  the  stability  boundary 
is  reached.  As  a  result,  the  methods  developed  in  this  work 
characterize  the  degree  of  confidence  of  the  computed 
stability  boundaries  by  examining  stability  under  finite 
external  disturbances.  A  first  order  memory  scheme,  which  can 
be  easily  implemented  in  real  time,  is  suggested  in  order  to 
expand  the  region  of  stability  of  straight  line  motion.  All 
computations  are  performed  for  the  NPS  autonomous  underwater 
vehicle  for  which  a  set  of  geometric  properties  and  slow 
motion  hydrodynamic  derivatives  is  available  [Ref.  1] .  Unless 
otherwise  mentioned,  all  results  are  presented  in  standard 
nondimens ional  form.  Equation  development  presented  in 
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[Ref. 13]  has  been  used  in  chapter  II,  and  with  modifications 
in  chapter  IV. 

C.  THESIS  OUTLINE 

Chapter  II  presents  the  formulation  of  vehicle  equations 
of  motion  and  the  criterion  for  determining  the  region  of 
guidance/control  stability  in  the  presence  of  a  position 
information  time  lag. 

Chapter  III  examines  the  vehicle's  control  stability 
through  the  use  of  a  Hopf  bifurcation  analysis. 

Chapter  IV  examines  the  effect  on  the  region  of  stability 
determination  by  the  use  of  a  memory  model  which  incorporates 
the  two  previous  vehicle  positions. 


II.  PROBLEM  FORMULATION 


A.  EQUATIONS  OP  MOTION 

The  mathematical  model  which  describes  vehicle  motion  in 
the  horizontal  plane  consists  of  nonlinear  sway  and  yaw 
equations  of  motion.  A  moving  coordinate  frame  which  is  fixed 
at  the  vehicle's  geometric  center  gives  Newtonian  equations 
of  motion  which  are 

m(  v+ui+xar-yji2)  =  Y  (2.1) 

Izi~mxG(v+ur)  -myGvz=N  (2.2) 


where  u  is  constant  forward  speed;  v  and  r  are  the  relative 
sway  and  yaw  velocities  of  the  moving  vehicle  relative  to  the 
water;  m  is  the  mass  of  the  moving  vehicle;  Xq,  yG  are  the 
respective  lateral  and  longitudinal  positions  of  the  center  of 
gravity;  and  Y,  N  represent  the  total  excitation  sway  force 
and  yaw  moment,  respectively.  The  vehicle's  added  mass, 
damping,  and  viscous  drag  due  to  motion  through  the  water  are 
represented  by  quadratic  drag  terms  and  memoryless  polynomials 
in  v  and  r.  Y  and  N  can  be  represented  as  the  sum  of  these 
terms  by 


y=JLj*Y t+JLl3  (y .vi-yur)  *^-l2Yvuv 
2*2  •'  *  2  ^ 


1=  ■£■  lsNti  +  -P- 1*  (N^v+NrUr)  +£-l3Nvuv 
2  *  2  **  *  2  *" 


where  1  is  the  vehicle  length,  Ya,  Na  represent  partial 
derivatives  of  Y  and  N  with  respect  to  a,  CDy  is  the  drag 
coefficient,  and  6  is  the  rudder  angle. 

At  relatively  high  speeds  and  low  angles  of  attack  the 
steering  response  is  predominantly  linear.  In  low  speed 
maneuvering  or  hovering  conditions  the  cross  flow  integral 
drag  term  becomes  significant. 


Vehicle  yaw  rate  is  expressed  by 


(2.3) 


and  inertial  position  rate  is  expressed  by 


y=  usim|f+vcosi|f 


(2.4) 


where  \p  is  the  heading  angle  and  shown  in  Figure  1. 
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Taking  equations  (2.1),  (2.2),  and  (2.3)  as  a  set  of  three 
coupled  linear  differential  equations,  and  a  linearized  form 
of  equation  (2.4),  the  final  set  of  equations  of  motion  for 
steering  control  are 


.  '* 


+ 


(2.5a) 


v,=  a11uv+a12ur+Jb1u28  (2.5b) 


f  =  a21uv+a22ur+b2u2t>  (2.5c) 


y =  uijr  + v  (2 . 5d) 

at  any  nominal  forward  speed. 

B.  GUIDANCE  AND  CONTROL 
1 .  Guidance  Scheme 

An  orientation  based  command  scheme  is  used  where  the 
commanded  vehicle  heading  angle  \J/C  is  a  function  of  the  line 
of  sight  angle  a  between  the  actual  vehicle  position  and  a 
reference  position  on  the  nominal  path  located  at  a  constant 
lookahead  distance  d.  Schematically,  this  is  presented  in 
Figure  1 . 

The  line  of  sight  angle  a  is  defined  by 
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tana =- — 
a 


(2.6) 


The  autopilot  will  be  called  upon  to  orient  the 
vehicle  to  the  commanded  heading  angle  \j/c,  which,  in  pursuit 
guidance,  will  equal  the  line  of  sight  angle.  This  defines  \pc 
as 

i|i  = -arctan-^  (2.7) 

a 


For  relatively  small  angles 


(2.8) 


2 .  Controller  Design 

The  steering  control  governing  equations  can  be 
represented  in  the  form 


x=Ax*bb  (2.9) 

where  the  state  vector  equation  is 

x-  (i|r,  v,  r]  T  (2.10) 


Linear  full  state  feedback  is  introduced  in  the  form 


6=kl(^~ifc)+k2v+kir  (2.11) 

The  gains  klt  k2,  and  k3  are  computed  by  pole  placement 
to  give  the  desired  dynamics  to  the  closed  loop  system.  The 
vehicle's  longitudinal  axis  is  pointed  toward  the  desired 
heading  by  the  difference  (\p-\pc)  in  the  control  law. 

With  a  dimensionless  time  constant  tc  the  general  form 
of  the  characteristic  equation  is 

A.3+a1Ai+a2l-a3  =  0  (2.12) 


where 


The  controller  gains  are  computed  from 

Jtx  = - ^ - 

(b2aii~bia2i)  u3 


(2.13a) 


^bXa22~b2ai2^  U  b2*  (b2ail  bia2\) 

=  az+bzu2kz- (aua„-al2a21)  u2 


(2.13b) 


bzu2kz*bzu2kz  =  -o1  -  (an+a22)  u 


(2.13c) 


These  gain  values  are  continuously  updated  due  to  changing 
nominal  forward  speeds.  Figure  2  shows  the  three  controller 
gains  at  unit  forward  speed  versus  the  system  time  constant. 

C.  TIME  LAG 

All  of  the  required  state  information  for  vehicle  control 
is  available  at  sufficient  rates  with  the  possible  exception 
of  the  y  position  information.  Due  to  time  requirements  for 
reduction  of  sonar  returns  and  inertial  navigation 
information,  there  may  be  a  time  lag  in  the  y  position 
information. 

This  time  lag  T  (sec)  is  introduced  to  the  guidance 
control  law  and  the  commanded  rudder  angle  \pc  becomes 


i|rc=-arctan>l^n 
c  d 

For  small  angles 

4,  s  -ziiiZL 

d 


(2.14) 


(2.15) 


♦ 


s«D 


+ 


•  • 


and  the  resulting  linearized  control  law  becomes 

8  =  XilzH  +k2v+k3r 


(2.16) 
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GAINS 


The  linearized  equations  of  motion  become 


'•*) 


+  =  r 


v  =  jb1u2/c1i^+  (a11u+b1u2ic2)  v 

♦  (a12u+i>jU2ic3)  r*b1u2-^y(  t-T) 


(2.17a) 


(2.17b) 


r  =  b2u2kiy+  ( a2lu+b2u2k2 )  v 

+  (a22u+b2u2k,)  r*b2u2—jy(  t-T) 


(2.17c) 


y  =  uijf  +  v 


( 2 . 17d) 


•  • 


D.  STABILITY  ANALYSIS 

Control  law  stability  is  determined  by  the  eigenvalues  of 
the  matrix  A  from  the  linear  system 

x- Ax  (2.18) 


with  the  state  vector 


*=  [♦,  v,z,y]T 


(2.19) 
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where  the  state  equations  have  been  linearized  around  straight 


line  motion.  The  characteristic  equation  of  (2.18)  is  a 
quartic  of  the  form 

k*+BV+Ck2+Dk+E=0  (2.20) 


where  the  coefficients  B,  C,  D,  and  E  are  functions  of  vehicle 
properties,  controller  gains,  and  lookahead  distance. 

The  characteristic  equation  will  give  a  pair  of  complex 
conjugate  roots  which  cross  the  imaginary  axis  under  the 
following  conditions. 

BCD-B2E-D2- 0  (2.21a)  • 


D 

B 


>  0 


(2.21b) 


Minimum  lookahead  distance  dcrit  for  stability  can  be 
determined  by  translating  these  conditions  to 

a1di+a2d+a, =0  (2.22a) 


d>  bia22-b2ai2~b2 
b2ail~bia21 


(2.22b) 


where 
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ai  S'ai«2-°3 


(2.23a) 


_  2#^)  ( ^1^22  ~^2ai2  ^2^ 

•^2ail“*fc>ia21 

- _ _ _a2u 

(b^-b.a^)  u 


(2.23b) 


-  (b1a22-b2al2-b2)  [Jb^-  (b1a22-jb2a12-Jb2)  u]  ct, 
(‘b2a11-bla21)  2u 


(2.23c) 


Analysis  of  the  system  when  operating  with  an  information 
time  lag  requires  approximation  by  either  a  first  order  Taylor 
expansion 


yU-T )  ±y-Ty 


(2.241 


a  second  order  Taylor  expansion 


y(t-T)  Ay-iy+llf 


(2.25) 


a  third  order  Taylor  expansion 


rp2  7*3 

y(t-T)  =  y~Ty+  ——y-  ——y 

Z  D 


(2.26) 


or  a  frequency  response  analysis  using  the  Nyquist  criterion. 
Figure  3  shows  the  region  at  stability  given  by  each  of  these 
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methods  with  a  one  second  (dimentionless)  time  lag.  It  can  be 
seen  that  the  agreement  is,  in  general,  very  good  among  the 
three  Taylor  series  approximations,  especially  as  the  time 
constant  tc  is  increased.  The  third  order  approximation 
coincides  with  the  exact  solution  from  the  frequency  response 
(Nyquist  criterion)  method.  The  first  order  method  requires 
the  highest  value  of  d  for  stability  at  a  given  tc,  and 
therefore  is  the  most  conservative  method  for  design  use.  For 
this  reason  the  first  order  method  is  used  in  the  next  chapter 
to  study  the  dynamics  of  the  system  as  the  critical  stability 
boundary  is  crossed. 
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Taylor  expansion  and  Nyquist  stability  analysis 
with  a  one  second  time  lag.  (Taylor  first 
order-solid,  Taylor  second  order-dashed,  Taylor 
third  order  and  Nyquist-dotted) 


III.  HOP?  BIFURCATION 


A.  INTRODUCTION 

It  can  be  numerically  confirmed  that  in  all  cases  of 
stability  loss  of  the  previous  chapter,  one  pair  of  complex 
conjugate  eigenvalues  of  the  corresponding  eigenvalue  problem 
crosses  transversely  the  imaginary  axis.  A  situation  like 
this  in  which  a  certain  parameter  is  varied  such  that  the  real 
part  of  one  pair  of  complex  conjugate  eigenvalues  of  the 
linearized  system  matrix  crosses  zero,  will  result  in  the 
system  leaving  its  steady  state  in  an  oscillatory  manner. 
This  loss  of  stability  is  called  Hopf  bifurcation  and 
generically  occurs  in  one  of  two  ways,  supercritical  or 
subcritical.  In  the  supercritical  case,  stable  limit  cycles 
are  generated  after  the  nominal  straight  line  motion  loses  its 
stability.  The  amplitudes  of  these  limit  cycles  are 
continuously  increasing  as  the  parameter  distance  from  its 
critical  value  is  increased.  For  small  values  of  this 
criticality  distance  the  resulting  limit  cycle  is  of  small 
amplitude  and  differs  little  from  the  initial  nominal  state. 
In  the  subcritical  case,  however,  unstable  limit  cycles  are 
generated  before  the  nominal  state  loses  its  stability. 
Therefore,  depending  on  the  initial  conditions  it  is  possible 
to  diverge  away  from  the  nominal  straight  line  path  and 
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converge  toward  a  different  steady  state  even  before  the 
nominal  motion  loses  its  stability.  In  many  cases  this  new 
steady  state  is  another  limit  cycle  of  considerably  higher 
amplitude.  This  means  that  in  the  subcritical  Hopf 
bifurcation  case  the  domain  of  attraction  of  the  nominal  state 
is  decreasing  and  in  fact  it  shrinks  to  zero  as  the  critical 
point  is  approached.  Random  external  disturbances  of 
sufficient  magnitude  can  throw  the  vehicle  off  to  an 
oscillatory  steady  state  even  though  the  nominal  state  may 
still  remain  stable.  After  the  nominal  state  becomes 
unstable,  a  discontinuous  increase  in  the  magnitude  of  motions 
is  observed  as  there  exists  no  simple  stable  nearby  attractors 
for  the  vehicle  to  converge  to.  Distinction  between  these  two 
qualitatively  different  types  of  bifurcation  is,  therefore, 
essential  in  design.  The  computational  procedure  requires 
higher  order  approximations  in  the  equations  of  motion  and  is 
the  subject  of  the  next  section. 

B.  COMPUTATIONS 

The  vehicle  equations  of  motion  are  written  in  the  form 

\|/  =  r  (3.1a) 

vr=anuv+a12ur+jb1u28  (3.1b) 


i  =  a21uv+a22ur  +b2  u  26 


(3.1c) 


y  =  usiruji  +  vcosijf 


;3.1d) 


In  order  to  capture  the  effect  of  rudder  saturation  we 
write  the  rudder  angle  6  in  the  form 


6  ~  ^ sat  tanh(yS-) 

°S«t 


where  6sat  is  the  saturation  limit  of  6,  typically  set  at  +0.4 
radians.  Equation  (3.2)  is  used  instead  of  a  hard  saturation 
function  because  of  its  smoothness,  which  is  important  for  the 
computation  in  this  section.  60  is  the  linear  control  law  of 
the  form 


80  =  icx  ( »|r  v+k,z 


(3.3) 


where 


tan\|rc 


_  _  y  ( t-T)  _  _  y-Ty 


using  the  first  order  approximation  for  y(t-T). 

The  state  equations  (3.1)  through  (3.4)  are  written  in  the 
compact  form 
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x-  f  (x) 


(3.5) 


where  f  (x)  is  a  nonlinear  function  of  the  state  variables 
vector 

x=  [i| i.v.r.y]  (3.6) 


Expanding  f  (x)  in  Taylor  series,  we  can  write  (3.5)  in  the 
form 


x  =  Ax+g(x)  (3.7) 

where  A  is  the  linear  system  matrix  and  g(x)  contains  all  the 
leading  nonlinear  terms  of  f (x) .  Due  to  the  port/starboard 
symmetry  of  our  problem,  g(x)  contains  only  third  order  terms. 
Defining  T  as  the  matrix  of  eigenvalues  of  A  and  applying  the 
transformation 

x  =  Tz  (3.8) 


'4 

-w 


4 


•  < 


the  state  equation  is  transformed  into  its  canonical  form 

£  =  T~1ATz+T~1g(  Tz)  (3.9) 


At  the  Hopf  bifurcation  point 
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(3.10) 


0 


T~lAT= 


K 


0 

0 


-w0  0  0 
0  0  0 
0  p  0 
0  0  g 


where  p  and  q  are  negative  and  u>0  is  positive.  In  the  new 
system  of  coordinates 


z  =  T~lX  (3.11) 

the  dynamics  of  the  system  are  generated  by  a  reduced  two 
dimensional  system  zt  and  z2,  since  the  coordinates  z3  and  z4 
correspond  to  the  eigenvalues  p  and  q  and  are  asymptotically 
stable.  These  stable  variables  z3,  z4  can  be  expressed  as  a 
function  of  the  critical  variables  z1(  z2;  these  functional 
expressions  are  at  least  of  third  order.  Therefore,  the 
stable  coordinates  z3,  z4  do  not  influence  the  third  order 
expansion  of  (3.7).  Using  the  above  expression,  we  can  write 
the  two  dimensional  system  in  z1#  z2  as 

=  -<ji0z2->-r1izl+r12z?z1+r11z1  z2*r14z2  (3.12a) 


=  wo  *i ♦  r2 1 . zl + . r22  zl  z2*r23  zx  z\  +  r2i  z\ 


(3.12b) 


where  the  coefficients  ri3  are  complicated  expressions  that 
are  derived  from  (3.9) 
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Equations  (3.12)  are  valid  at  exactly  the  Hopf 
bifurcation  point.  For  values  of  the  parameter  d  close  to  the 
bifurcation  point,  they  are  written  as 

=  a/ez1-  )  z2+F1  (zlt  z2)  (3.13a) 


i2  =  ((Oq+q'c)  z1+a'cz2+F2  (zltz2)  (3 . 13b) 

where  a  ',  o>  '  are  the  derivatives  of  the  real  and  imaginary 
parts  of  the  critical  pair  of  eigenvalues  evaluated  at  the 
bifurcation  point;  e  is  the  difference  of  the  parameter  d  from 
its  critical  value;  and  the  nonlinear  function  F2 ,  F2  remain 
the  same  as  in  (3.12), 

Fx{zr,z2)  =r11zi+z’12z?z2+r13z1z£+r14zj  (3.14a) 

F2(z1,z2)  =  z21zl  +r22z?  z2+z23z1z2  +r24z2  (3.14b) 

If  we  introduce  polar  coordinates  in  the  form 


zx  =i?COS0 


(3.15a) 


z2  =i?sin0 


(3.15b) 


equations  (3.13)  become 


R  =  a'eR+F1  ( R ,  0)  cos &+F2  (R,  0)  sin0 


(3 . 16a) 


rB-  (u>0+u>'t)  R+F2(R,B)cosB-F1(R,B)  sinB  (3.16b) 

Equation  (3.16a)  yields 

R  =  a'eR+P(B)  R3  (3.17) 

where  P  (0+2ir)  =  P  ( 6)  .  If  (3.17)  is  averaged  over  one  cycle  in 
6,  we  get  an  equation  with  constant  coefficients 

R  =  a'eR+KRl  (3.18) 

where 

K=—[2*P(B)dB  (3.19) 

2n  Jo 

Carrying  out  the  indicated  integration  in  (3.19)  we  obtain 

*■=-§  <3ri1+ru+r22+3r24)  (3.20) 

Existence  and  stability  of  limit  cycles  can  be  determined 
by  analyzing  the  equilibrium  points  of  the  averaged  equation 
(3.18),  which  correspond  to  periodic  solutions  in  z1#  z2  as 
can  be  seen  from  (3.15)  .  If  we  examine  equation  (3.18)  we  can 
see  that: 
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1)  If  a  '>0,  then: 

(a)  If  K>0  then  unstable  limit  cycles  coexist  with  the 
stable  equilibrium  for  £<0. 

(b)  If  K<0  then  stable  limit  cycles  coexist  with  the 
unstable  equilibrium  for  £>0. 

2)  If  a  '<0,  then: 

(a)  If  K>0  then  unstable  limit  cycles  coexist  with  the 
stable  equilibrium  for  £>0. 

(b)  If  K<0  then  stable  limit  cycles  coexist  with  the 
unstable  equilibrium  for  £<0. 

Therefore,  we  can  see  that  computation  of  the  above 
nonlinear  coefficient  K  can  distinguish  between  the  two 
different  types  of  Hopf  bifurcation: 

•  Supercritical  if  K<0 

•  Subcritical  if  K>0 
[Ref.  10] . 

C.  RESULTS 

Values  of  the  coefficient  K  have  been  calculated  for 
several  different  values  of  a  y  position  information  time  lag, 
over  a  span  of  system  time  constants,  using  the  Fortran 
program  presented  in  appendix  A.  This  program  has  also  been 
used  to  determine  the  vehicle  period  of  oscillation  in  the 
same  conditions.  Figure  4  shows  the  behavior  of  the  K 
coefficient  at  three  different  time  lags,  and  Figure  5  shows 
periods  of  oscillation  at  the  same  three  time  lags.  It  can  be 
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Figure  5.  Vehicle  oscillation  period  versus  tc  for 
values  of  time  lag  T.  (0.5  secasolid 
secadashed.  1.5  secadotted) 


seen  that  supercritical  bifurcations  are  ensured  for 
sufficiently  high  values  of  the  control  time  constant  tc.  For 
tc  less  than  a  certain  critical  value,  the  corresponding  Hopf 
bifurcations  are  subcritical.  This  situation  should  be 
avoided  in  practice. 

D.  SIMULATIONS 

The  vehicle  path  has  been  simulated  by  using  the  Euler 
integration  method  coded  in  a  Fortran  program  presented  in 
appendix  B.  The  vehicle  control  law  (2.17)  as  well  as  the 
control  gains  (2.13)  are  used.  The  program  has  been  run  with 
input  values  for  system  time  constant,  y  position  information 
time  lag,  and  lookahead  distance.  The  vehicle  was  given  an 
initial  nominal  y  offset  of  half  a  shiplength  and  a  nominal 
forward  speed.  The  resulting  plots  of  y  position  against  time 
show  the  vehicle’s  oscillatory  path  either  converge  to  the 
commanded  path  or  diverge. 

The  results  of  these  simulation  runs  were  compared  to  the 
results  of  the  Hopf  bifurcation  computations  in  two  ways.  The 
period  of  oscillation  observed  in  the  vehicle  path  was 
compared  to  the  period  predicted  by  the  Hopf  program. 
Additionally,  the  vehicle  stability,  determined  by  convergence 
to  the  commanded  path,  was  compared  to  the  K  coefficient  sign 
prediction  of  stability  which  was  determined  by  the  Hopf 
computations.  Figure  6  shows  a  vehicle  path  in  stable 
conditions.  Figure  7  shows  a  vehicle  path  in  unstable 


28 


conditions.  These  simulations  were  chosen  in  the 
supercritical  bifurcation  case  and  it  can  be  seen  that  the 
period  of  oscillation  observed  from  Figures  7  agrees  very  well 
with  the  theoretical  results  of  Figure  5. 


Figure  6.  Vehicle  path  in  stable  conditions  where  time 
constant  tc  ■  1.0  sec,  time  lag  T  ■  1.0  sec,  and 
lookahead  distance  d  -  2  shiplengths. 


IV.  STABILITY  ENHANCEMENT 

A.  TWO- POINT  FORMULA 

The  time  lag  of  y  position  information  discussed  in 
chapter  II  has  thus  far  been  represented  by  a  single  piece  of 
information  used  in  the  control  law  with  some  delay  assigned 
to  it.  In  an  attempt  to  improve  stability  with  regard  to  this 
delay,  the  use  of  the  previous  two  positions  in  the  control 
was  examined. 

Equation  (2.15)  shows  the  representation  of  the  time 
delay,  and  equation  (2.16)  shows  its  effect  on  the  control 
law.  A  straight  line  approximation  is  applied  to  the  previous 
two  positions  as  shown  in  Figure  8.  The  general  y  position  is 
defined  by 

y  ( t)  =  a1c+a2  (4.1) 

The  two  previous  positions  become 

y1=a1t1+a2  (4.2a) 


* 


£ 
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y2  -  a1 t2  +a2 

The  coefficients  ax  and  a2  are  stated  as 


(4.2b) 
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*  -  yi'Y* 

d,  —  * — — - 

1  tx-c2 


(4.3a) 


a2 


y^z-yzti 

t2-ti 


(4.3b) 


The  previous  time  values  are  represented  as 

t.  *  C-2T  (4.4a) 


t2  =  t-T 


(4.4b) 


where  T  is  the  time  lag  associated  with  the  y  position  ^ 

information.  The  previous  positions,  in  terms  of  the  previous 
time  values  are 

yx=y(C-2T)  (4.5a)  • 


y2=y(t-T) 


(4.5b) 


Substitution  and  algebra  gives  a  general  position 
expression  of 


y  ( t)  =2  y(t-T)  -y(t-2T) 


(4.6) 


•  •  •  •  • 
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B.  STABILITY  ANALYSIS 


The  Laplace  transform  of  equation  (4.6)  is 

y(t)  =y  (2e'Ts-e'2Ts)  (4.7) 


The  control  law  becomes 

6  =  *k2v+k3r  +  ~ y  (2e~Ts-e~2Ts) 


(4.8) 


and  the  linearized  equations  of  motion  become 


ij,  =  r 


(4.9a) 


v  =  jb1u2ic1i|r  +  (a11u+b1u2k2)  v+  (a12u+bxu2k2)  r 
+b1u2  ■^ky(2e~Ts-e~2Ts) 


(4.9b) 


r  -  £>2u2/c1i|r  +  ( a12u+b2u2k2 )  v+  (a22u+b2u2k2)  r 

kx 
~d 


+£>,  u 2  ^1  y  ( 2  e  "Ts- e  ~2Ts) 


(4.9 c) 


y  =  ui|r  +  v 


(4 . 9d) 


These  motion  equations  reduce  to  the  state  space  form 

x-Ax  (4.10) 
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The  matrix  form  becomes 


v  b3u2kx  a11u+jb1u2ic2  a12u+b3u2k3  bxu2^±  (2e'rs-e'2Ts)|v 
i  b2u2k3  a21u+b2u2k2  a22u+b2u2k3  b2u2^  (2e~Ts-e~2Ts)\r 


The  characteristic  equation  of  the  form  [A-Is]=0  is 


s4+A3s3+A2s2+A1s+  (B2s2*bxs*B0)  D(2e'Ts-e-2Ts)  (4.11) 


where 


A,  =  -  (a11+a22)  v.  -  ( b1k2*b2k3 )  u 2 


A2=  (*11*22 -*12*21)  u2*(a31b2-a21b3)  u3ic3  +  (a22b3-a12b2)  u3k2-b2u'-k3 


Aix~(a21b1-allb2)  u3iq 


B2  =  -blu2k1 


B1  =  ~b2 u2k3-{ a12b2 -a22bx )  u3k3 


•  4 


I  ( mmk 


B0=(aiA-a2A)  U'kl 


TT 


'M 
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The  characteristic  equation  of  the  system  is  written  as 

l+*G(s)  =0  (4.12) 


where 


G(S)  = 


(B2s2+Bxs+B0)  {2e'Ts~e'2Ts) 
S  ( s3  -rAjS2  ■>-A2s+A1) 


(4.13) 


is  the  transfer  function,  and  we  have  denoted 


•  • 


K=  D 


(4.14) 


With  s=joj,  the  phase  angle  is  represented  by 

4>  =  ZG(J0))  (4.15) 


where 


<j>  =  i  (2e'rju-e'2T;“)  -Z  ( jo))  +Z  (  -B2(j>2+B1ju>  +B0) 
-L  (-ju>3~A3u>2*A2j(t>*Al) 
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4>  =  tan'1  -2sin((i>r)  -t-sin  (_2<j_T)  _*  +tan-i  ( 

2cos  (uD  -cos  (2o»D  2  B0-B2u> “ 

.  -u3-A2u>  . 

-tan  1  ( - r-^— ) 


(4.16) 


The  Nyquist  stability  criterion  states  that  at  the  solution  of 


4>  (g>)  =  -it 


(4.17) 


where  u  is  the  phase  cross  over  frequency  u)1(  the  gain  margin 


must  equal  1 


| KGljiaJ  |=1 


where  K  =  D  =  1/d.  Equation  (4.17)  becomes 


(4.18) 


•  • 


-2sin (o^T)  +sin Uci^T)  n 


2cos(g)1D -cos  2 


i-tan'M- 


b0-b2oi 


.  "Uj+AjWj  . 

-tan  1  ( — -  )  =  -n 


Rearranged  as 


(4.19) 


.  .  B.o).  %  .  ,  -Oi+A,0).  4 

tan  1  ( - i— i— )  -tan'1  ( -  2  1 )  = 

B0~B2^i 


n  ,  -2sin(o),D +sin(2(i),T) 

-  —  -  tan  1 - - - - — 

2  2cos  ((02T)  -cos  (2o>1T) 


Taking  the  tangent  of  both  sides  and  using  the  identity 


f: 


tan (x-y) 


tan(x) -tan (y) 
l-tan(x) tan(y) 


then  algebra  and  rearranging  gives 


Bi<*i  (Ax-A3<i>i)  -  ( B0-B2(jil )  (-<^1*A2<j)1) 

(B0-B2u>\)  (-A3<*1+Ax)  +S1q1(-Ui+A2w1)  (4.20) 

2cos  (u1T-cos  (2 u>xT) 

-2sin(<d1r)  +sin(2w1D 


£ 


The  stability  computations  require  an  initial  value  for 
This  is  accomplished  by  setting  the  time  lag  T=0,  and 

using 

(Bq-B2u>\)  i-A^l+A^  *B1oi1{-Uil+A2(ti1)  =0  (4.21) 


•  • 


rearranged  as 


( B2A-,  -B1 )  co  1  +  ( BXAZ -B2A 1  -B0A3 )  <j>\+B0Ax  =  0 


(4.22) 


Setting  the  gain  margin  equal  to  1  gives 
1  |(-BaM;+B1jM1-fB0)  (2e~j'*'T-e'2:i'*'T)\ 


and  d  can  be  solved  for  by 
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(Ai~A3<*1)  2  +  (A2a>1-<iil) 2 


(4.23) 


In  order  to  facilitate  computation,  terms  are  grouped  as 
follows 


m 


4r 


*1 


«z  =S1A2-S2A1-B0Aj 


a,  =B0A. 


Ix—B2 


P2=B0*fi2A2-fiiAi 


P3=b1a-b0^2 


Equation  (4.20)  becomes 

P1fa>5^P2to3»P3fa>  _  2cos  (t)D  -cos  (2ti)D 
o1(«)4+a2ti)2+o3  -2sin(wT)  -*-sin(2«D 


» 


Put  in  the  form 
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equation  (4.24)  becomes 


[  (P1w5>P2uJ'*-P3«)  (-2sin(<i)D  +sin(2 uD  )  ] 

-  [  (a1<i)4>a2a>2-*-a3)  (2cos  (ul)  -cos  (2 uD  )  ]  =0 


Newton's  iteration  is  used  as 


(»*-!> 


Where  the  function  of  w  is 

f(u)  =  01<i>5  +  P2a>3+P3<i>)  (-2sin(wT)  ) 
♦  (P1cjs+P2w3-t-P3o)>  (sin(2c«>r) ) 
-  (a1<j4+«2w2+o3)  (2cos(o>D) 

+  (a1cj4+a2«2+o3)  (cos  (2  wD) 


and  the  derivative  function  is 

/'(«)  =  (5p.w4+3P2<i)z+P3)  ( -2sin («D  ) 

-  (  P10)s  +  P2u)3+P1u)  (2 TcositoT) 

+  (5P1oj*+3P2o)2  +  P1)  (sin  (2«D  ) 

+  (P1w5+P2o>3-»-P3(i))  (2Tcos(2fc)D) 
+  (4a1w3+2o2w)  (~2cos(uD) 

+  («1Ci)4+«2o)2+o3)  (2rsin(wT)) 

+  (4o1«3+2o2«)  (cos  (2«T)  ) 

-  (o1ci)4+a2w2+o3)  (2Tsin(2wT) ) 


(4.26) 


(4.27) 


(4.28) 


(4.29) 


Appendix  C  presents  the  Fortran  program  used  to  determine 
stability  characteristics  as  a  function  of  system  time 
constant  and  time  lag. 

C.  RESULTS  AMD  DISCUSSION 

Values  for  the  minimum  lookahead  distance  required  for 
controller  stability  have  been  computed  as  a  function  of  the 
system  time  constant  and  the  y  position  information  time  lag 
using  the  two-point  formula.  Figures  9,  10,  and  11  show 
critical  lookahead  distance  plotted  against  time  constant  at 
three  values  of  time  lag  for  both  the  two -point  method  as  well 
as  the  single  point  (Nyquist)  method  discussed  in  chapter  II. 
These  figures  show  an  overall  decrease  in  required  lookahead 
distance  as  a  result  of  the  use  of  the  two-point  memory  model. 
An  exception  to  this  occurs  at  low  values  of  tc.  However, 
these  values  are  to  be  avoided  in  practice,  as  explained  in 
the  previous  chapter. 

Figure  12  through  15  show  the  critical  lookahead  distance 
plotted  against  time  lag  for  different  values  of  the  time 
constant  using  both  the  two-point  and  single  point  methods. 
Again  the  two-point  method  produces  superior  results  to  the 
single  point  method,  with  the  exception  of  small  time 
constants  and  large  time  lags.  The  same  information  is 
summarized  in  Figure  16  through  19,  where  we  show  the  ratio 
d2/dx  (critical  lookahead  distance  using  the  two-point  method 
divided  by  the  critical  lookahead  distance  using  the  single 


point  method)  verses  time  constant  and  time  lag.  It  can  be 
seen  that  the  use  of  the  two -point  method  can  result  in 
improvement  in  the  region  of  stability  of  over  20%. 


0.8  1  1.2  1.4  1.6  1.8  2  2.2 

d 

Figure  9.  Critical  lookahead  distance  d  versus  time  constant 
tc  with  time  lag  T  »  0.5  sec  (two-point 

method>solid,  single  point»dashed) . 
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1.6 


d 

Figure  11.  Critical  lookahead  distance  d  versus  time  constant 
tc  with  time  lag  T  «  1.5  sec  (two-point 

method-solid,  single  point-dashed) . 
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0.75  0.0  0.85  0.9  0.95  1  1.05  1.1  1.15  l.Si 

d2/dl 


Figure  16.  Ratio  dj/d^  versus  time  constant  tc  with  time 
lag  ■  1.0  sec. 


d2/dl 


Figure  17.  Ratio  d2/dx  versus  time  constant  tc  with  time 
lag  -  1.5  sec. 


figure  18.  Ratio  d^d-^  versus  tine  lag  T  with  time  constant  » 
1.0  sec 
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V. 


CONCLUSIONS  AND  RECOMKKNDAT IONS 


A.  CONCLUSIONS 

Nonlinear  bifurcation  theory  has  been  applied  to  the  NPS 
autonomous  underwater  vehicle  through  a  Hopf  bifurcation 
analysis  of  the  vehicle's  navigation  guidance  and  control 
scheme.  It  has  been  shown  that  under  normal  conditions  the 
vehicle  will  operate  in  the  supercritical  region,  and 
therefore  the  unpredictable  behavior  associated  with  the 
subcritical  case  should  not  to  be  expected. 

Secondly,  the  control  scheme  was  modified  by  the  use  of  a 
single  stage  memory  model  which  incorporates  the  two  previous 
vehicle  positions  in  the  linearized  control  law.  It  has  been 
shown  that  the  use  of  this  memory  model  reduces  the  vehicle's 
required  lookhead  distance  for  control  stability  within  the 
normal  operating  range. 

B.  RECOMMENDATIONS 

In  the  event  that  not  all  of  the  state  information 
required  is  directly  available  through  measurements,  the 
control/guidance  scheme  would  include  an  observer.  The  effect 
of  the  observer  dynamics  on  the  higher  order  nonlinear  terms, 
and  therefore  on  the  bifurcation  limit  cycle  stability,  should 


be  examined. 


APPENDIX  A 


C 

C 

C 

C 

C 

C 

C 


C 


C 


C 

C 

C 


C 

C 

C 


PROGRAM  HOPF. FOR 

Hopf  Bifurcation  Analysis 

Critical  Points  Exact 

Third  Order  Expansions:  rirst  Order  Approximation 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DOUBLE  PRECISION  Kl , K2 , K3 , K4 , L , NR , NV , NDRS , NDRB , I 2 , MASS , 

*  NRDOT,NVDOT,KlP,K2P 

DOUBLE  PRECISION  Mil , M12 , Ml 3 , M14 , M21 , M22 , M2 3 , M24 , 

1  M31,M32,M33,M34,M41,M42,M43,M44, 

2  N11,N12,N13,N14,N21,N22,N23,N24, 

3  N31,N32,N33,N34,N41,N42,N43,N44, 

4  L21,L22,L23,L24, L31 , L32 , L33 , L34 , 

5  L41,L42,L43,L44 

DIMENSION  A( 4 , 4 ) , T( 4 , 4 ) , TINV( 4 , 4 ) , FVl ( 4 ) , IVl ( 4 ) , ZZZ ( 4 , 4 ) 
DIMENSION  WR( 4 ) ,WI ( 4 ) ,TSAVE( 4,4), TLUD<  4,4), IVLUD( 4 ) , SVLUD( 4 ) 
DIMENSION  ASAVE (4,4) 

OPEN  (11, FILE- ' HOPF . DAT ' , STATUS- ' NEW ' ) 

OPEN  (12,FILE-'PERI .DAT' , STATUS- ' NEW ' ) 

OPEN  (13, FILE- 'ALPH10.DAT' , STATUS- ' NEW' ) 

Vehicle  Parameters 

WEIGHT-435.0 
IZ  -45.0 

L  *7 . 3 

RHO  -1.94 

G  -32.2 

XG  -0.0104 
MASS  -WEIGHT/G 
U  -5.0 

RATIO  -0.0 

DO  -0.4 

WRITE  (*,1006) 

READ  ( * , *  )  DO 

YRDOT  — 0.00178*0. 5*RHO*L**4 
YVDOT  — 0.03430*0. 5*RHO*L**3 
YR  -+0. 01187*0. 5*RHO*L**3 
YV  — 0.03896*0. 5*RHO*L**2 
YDRS  -+0. 02345*0. 5*RHO*L**2 
YDRB  -+0. 02345*0. 5*RHO*L**2 
NRDOT  — 0.00047*0. 5*RHO*L**5 
NVDOT  — 0.00178*0. 5*RHO*L**4 
NR  — 0.01022*0. 5*RHO*L**4 
NV  — 0.00769*0. 5*RHO*L**3 
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NDftS  —  0.337*0'. 02345*0. 5*RHO*L**3 
NDRB  -+0.283*0.02345*0. 5*RHO*L** 3 


best 


OH  — ( IZ-NRDOT )  *  ( MASS- YVDOT ) - 
i  ( MASS*XG-YRDOT ) * ( MASS  *XG-NVDOT ) 

AAl 1— ( ( I Z-NRDOT ) *  YV- ( MASS*XG-YRDOT ) *NV ) /DH 
AA12- (  ( IZ-NRDOT)  M-MASS+YR)- 
i  (MASS*XG-YRDOT)*(-MASS*XG+NR) )/DH 

AA21- ( ( MASS-YVDOT ) *NV- ( MASS*XG-NVDOT ) * YV ) /DH 
AA22- ( ( MASS-YVDOT ) * ( -MASS*XG+NR ) - 
i  ( MASS*XG-NVDOT ) * ( -MASS+YR ) )/DH 

BB11- ( ( IZ-NRDOT ) *YDRS- ( MASS *XG- YRDOT ) *NDRS ) /DH 
BB12- ( ( IZ-NRDOT) *YDRB-(MASS*XG- YRDOT) *NDRB)/DH 
BB21- ( ( MASS-YVDOT ) *NDRS- ( MASS*XG-NVDOT ) * YDRS ) /DH 
BB22- ( ( MASS-YVDOT ) *NDRB- ( MASS*XG-NVDOT ) *  YDRB ) /DH 


Available  Cop) 


BBl— BBll+RATIO*BBl2 
BB2— BB21+RATIO*BB22 

WRITE  (*,1005) 

READ  ( * , * )  TL 

WRITE  (*,1001) 

READ  ( * , *  )  TC 


EPS-1. D-10 
ITMAX-1000 

TL-0.5 

DO  50  M-1,100 
TLD-TL*L/U 

TC-0.5 

DO  40  K-1,100 
TCD-TC*L/U 
POLEC— 1 . 0/TCD 
ADI-3 . 0*POLEC 
AD2-3 . 0* POLEC* *2 
AD3-POLEC**3 
Al-BBl *U*U 
B1-BB2*U*U 

Cl— AD1-(AA11+AA22  )*U 

A2- ( BBl *AA22-BB2*AAl 2 ) *U*  *  3 

B2*( BB2*AAll-BBl*AA21 ) *U**3 

Kl-AD3/( ( BB2  *AAl 1-BBl *AA21 ) *U*  *  3 ) 

C2-AD2-(AA11*AA22-AA12*AA21 ) *U**2+BB2*U*U*Kl 

K2-(C1*B2-C2*B1 )/(Al*B2-A2*Bl ) 

K3-(C2*A1-C1*A2 )/( Al*B2-A2*Bl ) 

Al-  ( AAl 1 *BB2-AA21 *BBl )*U*U*U*Kl 

A2-  (AA11*AA22-AA12*AA21 ) *U*U+ ( AAl 1 *BB2-AA21 *BBl ) *U*U*U*K3+ 
&  (AA22*BB1-AA12*BB2 ) *U*U*U*K2-BB2*U*U*Kl 

A3— (AA11+AA22 )*U-(BB1*K2  +  BB2*K3  )  *U*U 
B0-  (AA11*BB2-AA21*BB1 )*U*U*U*U*K1 
Bl— (AA12*BB2-AA22*BB1 )*U*U*U*K1-BB2*U*U*U*K1 
B2— BBl*U*U*Kl 

Compute  "oitial  Approximation  For  Omega  (TL-0) 


AQ-B2*A3-B1 

BQ-Bl*A2-82*Al-B0*A3 

CQ-B0*A1 

DET-BQ**2-4.D0*AQ*CQ 

IP  (DBT.LT.O.DO)  GO  TO  40 

ZTl-  ( -BQ+DSQRT { DET ) )/(2.0D0*AQ) 

ZT2- ( -BQ-DSQRT( DET ) )/(2.0D0*AQ) 

IP  (ZT1.GE.0.D0)  OM-DSQRT(ZTl) 

IP  (ZT2.GE.0.D0)  OM-DSQRT(ZT2 ) 

ALPHA1-AQ 
ALPHA2-BQ 
ALPHA3-CQ 
BETA1— B2 

BETA2»B0+B2*A2-Bl*A3 

BETA3-B1*A1-B0*A2 

Compute  Exact  Omega  Using  Newton's  Iteration 
OMOLD-OM 

P-(BETA1*OMOLD**5+BETA2*OMOLD**3+BETA3*OMOLD)*DSIN(OMOLD*TLD) 
+(ALPHAl*OMOLD**4+ALPHA2*OMOLD**2+ALPHA3)*DCOS(OMOLD*TLD) 
PPRIME-( 5.D0*BETAl*OMOLD**4+3 .DO*BETA2*OMOLD**2+BETA3 ) 
*DSIN(OMOLD*TLD)+(BETAl*OMOLD**5+BETA2*OMOLD**3+BETA3*OMOLD) 
*TLD*DCOS ( OMOLD+TLD )  +  (  4  . D0+ALPHA1 *OMOLD*  *  3  +  2 . DO  *ALPHA2*OMOLD ) 
+  DCOS ( OMOLD*TLD ) - ( ALPHAl +OMOLD*  *  4+ALPHA2  *OMOLD*  *  2+ALPHA3 ) 
*TLD*DSIN(OMOLD*TLD) 

IF  (FPRIME.EQ.0.D0)  STOP  1112 
IF  (F.EQ.0.D0)  GO  TO  92 
OMNEW-OMOLD-F/FPRIME 
OMDI F-DABS ( OMNEW-OMOLD ) 

IF  (OMDIF.LE.EPS)  GO  TO  92 

OMOLD-OMNEW 

IT-IT+1 

IF  ( IT.GT. ITWAX)  STOP  1111 
GO  TO  93 

OM-OMNEW 

XDNUM-DSQRT( ( B0-B2 *OM* * 2 ) *  * 2+ ( B1 *OM ) * * 2 ) 

XDDEN«OM*DSQRT( (Al-A3*OM**2 ) **2+< A2*OM-OM**3 ) **2 ) 
XD-XDNUM/XDDEN 

Start  Hopf  Bifurcation  Analysis 

Evaluate  Nonlinear  Rudder  Expansion  Coefficients 

K1P-K1-!U*TLD*U/XD 

K2P-K2-K1*TLD/XD 

DPPV—  (l.D0/(3.D0*D0**2) )*3.D0*K1P*K1P*K2P 

+  0 . 5D0*Kl*TLD/XD  +  3.D0*Kl*(TLD*U)**3/( 3.D0*XD**3) 

DPW—  (l.D0/(  3.D0*D0**2)  )  *3  .  D0*K1P*K2P*K2P 
+  3 . DO*U*TLD*TLD*TLD*Kl/( 3 . D0*XD*  *  3 ) 

DPPR—  ( 1 . D0/( 3.D0*D0**2) ) *3 .D0*K1P*K1P*K3 
DPRR—  ( 1 .  D0/(  3  .  DO* DO* *2  )  )  *3 .  D0*K1P*K3 *K3 
DPPY—  ( 1 . D0/( 3.D0*D0**2) ) *3 .D0*K1P*K1P*K1/XD 
-  3 . D0*TLD*TLD*U*U*Kl/( 3 . D0*XD** 3 ) 

DPYY—  ( 1 .D0/( 3.D0*D0**2 ) ) *3 . D0*K1P*K1 *Kl/( XD* *2 ) 

+  3 . D0*TLD*U*Kl/( 3 . D0*XD*  *  3 ) 

DWR—  ( 1  .D0/<  3.D0*D0**2)  )  *3  .  D0*K2P*K2P*K3 
DVRR—  (l.D0/(3.D0*D0**2) ) *3 . D0*K2P*K3*K3 


DWY—  ( 1 .  D0/(  3  .  D0*D0**2  )  )*3.D0*K1*K2P*K2P/XD 
t  -  3 . DO*Kl*TLD*TLD/( 3 . DO*XD*  *  3 ) 

DVYY—  ( 1 .D0/{ 3.D0*D0**2) ) *  3 . DO*Kl *Kl*K2P/( XD* *2 ) 

&  +  3 . DO*TLD*Kl/( 3 . DO*XD*  *  3 ) 

DRRY—  ( 1 ,D0/( 3 .D0*D0**2 ) ) *3 .D0*K1*K3*K3/XD 
DRYY—  ( 1 .  D0/( 3 .DO* DO* *2 ) )*3.D0*Kl*Kl*K3/( XD*  *2 ) 

DPVR—  ( 1 .D0/( 3.D0*D0**2 ) ) *6 . D0*K1P*K2P*K3 
DPVY—  ( 1 .D0/( 3.D0*D0**2  )  )  *6  .  DO ‘KIP’*  K1  *K2P/XD 
&  -  6.D0*TLD*TLD*U*Kl/( 3.D0*XD**3) 

DPRY—  ( 1 .  DO/(  3 . DO* DO*  *2 )  )  *6  .  DO*KlP*Kl  *K3/XD 
DVRY—  (l.D0/(3.D0*D0**2) ) *6 .D0*K1*K2P*K3/XD 
DPPP—  ( 1 . DO/( 3.D0*D0**2) )*1.D0*K1P*K1P*K1P 
t  +  Kl*TLD*U/( 6 . D0*XD )  +  ( Kl* ( TLD*U ) *  *  3 )/ ( 3 . D0*XD* *  3 ) 

DWV—  <l.D0/(  3 . D0*D0**2 )  )  *1 .  D0*K2P*K2P*K2P 
«  +  Kl * ( TLD**  3 )/( 3 . DO  *XD*  *  3 ) 

DRRR—  (l.D0/(3.D0*D0**2) )*1.D0*K3*K3*K3 

DYYY—  ( 1 .D0/( 3.D0*D0**2 ) ) * ( Kl/XD) **3-Kl/( 3 .D0*XD**3 ) 

&  -  Kl/( 3 . D0*XD* * 3 ) 

Evaluate  Transformation  Matrix  of  Eigenvectors 

PSI  -0.D0 

Y  -0.D0 

V  -0.D0 
DPSI-K1P 
DV  -K2P 
DR  -K  3 

DY  -Kl*XD/(XD**2+Y**2) 

A(1,1)-0.0D0 

A(1,2)-0.0D0 

A(1,3)-1.0D0 

A(1,4)-0.0D0 

A( 2 , 1 )  -  BB1*U*U*DPSI 

A(2,2)-AA11*U+BB1*U*U*DV 
A(  2 , 3 )-AA12*U+BBl*U*U*DR 
A( 2 , 4  )-  BBl *U*U*DY 

A( 3 , 1  )-  BB2*U*U*DPSI 

A( 3,2)-AA21*U+BB2*U*U*DV 
A(  3 , 3 )-AA22*U+BB2*U*U*DR 
A(3,4)-  BB2*U*U*DY 

A( 4, 1 )-U*DCOS( PSI )-V*DSIN( PSI  ) 

A( 4 , 2  )-  DCOS(PSI) 

A(4,3)-0.0D0 
A(4,4)-0.0D0 
DO  11  1-1,4 
DO  12  J-1,4 

ASAVE ( I ,  J  ) -A ( I , J ) 

12  CONTINUE 

11  CONTINUE 

CALL  RG( 4 , 4 ,A,WR,WI , 1 ,  ZZZ , IVl , FVl , I  ERR ) 

CALL  DSOMEGi IEV ,WR ,WI , OMEGA , CHECK ) 

OMEGAO-OMEGA 
DO  5  1-1,4 

T(I,1)-  ZZZ ( I , I EV ) 

T(  I  ,  2  )— ZZZ(  I  ,  IEV+1 ) 

5  CONTINUE 

IF  (IEV.EQ.l)  GO  TO  13 
IF  (IEV.EQ.2)  GO  TO  14 
IF  (IEV.EQ.3)  GO  TO  15 
STOP  3004 
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DO  6  1-1,4 

T(  1 , 3 )— ZZZ (1,1) 

T(  1 , 4 )-ZZZ (1,4) 

CONTINUE 
GO  TO  17 
DO  7  1-1,4 

T(  1 , 3 )-ZZZ (1,1) 

T(  1 , 4  )-ZZZ (1,2) 

CONTINUE 
GO  TO  17 
DO  16  1-1,4 

T( I , 3 )-ZZZ (1,3) 

T( I , 4 )-ZZZ (1,4) 

CONTINUE 

CONTINUE 

Noraalization  of  the  Critical  Eigenvector 

CALL  NORMAL ( T ) 

Definition  of  Mij 

Ml 1— T (1,1) 

M21-T( 2,1) 

M31«T( 3,1) 

M41-T( 4,1) 

M12-T( 1,2) 

M22-T( 2,2) 

M32-T( 3,2) 

M42-T( 4,2) 

M13«*T(  1,3) 

M23-T( 2,3) 

M33-T( 3,3) 

M43-T<  4,3) 

M14-T( 1,4) 

M24-T( 2,4) 

M34-T( 3,4) 

M44-T<  4,4 ) 


Definition  of  Lij 

L21-  DPPV*Mll*Mll*M21  +  DPVV*M11*M21*M21 
+  DP*R*Mll*M31*M31  +  DPPY*Mll*Mll*M41 
+  DWR*M21*M21*M31  +  DVRR*M2 1  *M3 1  *M3 1 
+  DVYY*M21*M41*M41  +  DRRY*M3 1 *M3 1 *M4 1 
+  DPVR  *M11*M21*M31  +  DPVY*Mll*M21*M41 
+  DVR Y  *M21*M41*M31  +  DPPP*Mll*Mll*Mll 
+  DRRR*M31*M31*M31  +  DYYY*M41*M41*M41 


DP PR* Ml l*Mll*M31 
DPYY*M11*M41*M41 
DVVY*M2 1 *M2 1 *M4 1 
DRYY*M31 *M4 1 *M4 1 
DPRY*M11*M41*M31 
DVW*M21*M21*M21 


Ml 1 *Ml 1 
M12*M21 
Mil *Ml 1 
M12*M31 
Mll*Mll 
M41*M41 
M21 *M21 
M22*M31 
M21 *M21 
M22*M4 1 
M31*M31 


*M22  + 
*M21  + 
*M32  + 
*M31  + 
*M42  + 
*Ml2  + 
*M32  + 
*M31  + 
*M42  + 
*M41  + 
*M42  + 


2 . 0*M11 
2 . 0  *Ml 1 
2 . 0*M11 
2 . 0*M1 1 
2 . 0*Mll 
2 . 0*M1 1 
2 . 0*M31 
2 . 0*M31 
2 . 0*M41 
2 . 0*M41 
2 . 0*M41 


*Ml2*M21 
*M21*M22 
*M12*M31 
*M31*M32 
*M12*M41 
*M4 1 *M4  2 
*M21*M22 
*M32*M21 
*M21*M22 
*M42*M21 
*M31*M32 


Best  Available  Cop* 


RYY  -  M32*M41*M41  +  2 . 0*M41*M42*M31 
PVR  -  Mll*M21*M32+M31*(Mll*M22+Ml2*M21 ) 

PVY  -  Mll*M21*M42+M41MMll*M22+Ml2*M21 ) 

PRY  -  Mll*M41*M32+M31*(Mll*M42+Ml2*M41 ) 

VRY  -  M21*M41*M32+M31*(M21*M42+M22*M41) 

PPP  -  3.0*M11*M11*M12 
VW  -  3.0*M21*M21*M22 
RRR  -  3.0*M31*M31*M32 
YYY  -  3.0*M41*M41*M42 

L22-DPPV*PPV+DPW*PW+DPPR*PPR+DPRR*PRR  +  DPPY*PPY+DPYY*PYY 
fc  +  DVVR*  WR+DVRR*VRR+DWY*  WY+DVYY*  VYY  +  DRRY*RRY+DRYY*RYY 

&  +DPVR*PVR+DPVY*PVY+DPRY*PRY+DVRY*VRY+DPPP*PPP+DVVV*VVV 

fc  +DRRR*RRR+DYYY* YYY 


M12*M12 
Ml 1 *M22 
M12*M12 
Ml 1 *M32 
M12*M12 
Mil *M42 
M22*M22 
M21 *M32 
M22*M22 
M21*M42 
M32*M32 
M31*M42 
M12*M22 
M12*M22 
M12*M42 
M22*M42 
3 . 0*M1 1 
3.0*M21 
3 . 0*M31 
3 . 0  *M4 1 


*M21 

*M22 

*M31 

*M32 

*M4 1 

*M42 

*M31 

*M32 

*M41 

*M4  2 

*M41 

*M42 

*M31 

*M41 

*M31 

*M31 

*M12* 

*M22* 

*M32* 

*M42* 


+  2.0 
+  2.0 
+  2.0 
+  2.0 
+  2.0 
+  2.0 
+  2.0 
+  2.0 
+  2.0 
+  2.0 
+  2.0 
+  2.0 
+  M32 
+  M42 
+  M32 
+  M32 
Ml  2 
M22 
M32 
M42 


*Ml  1  * 
*M12* 
*Ml  1  * 
*M12* 
*Ml  1  * 
*M12* 
*M21* 
*M22* 
*M2 1  * 
*M22* 
*M31* 
*M32* 
M  Mil 

*  ( Ml  1 
MM11 

*  ( M21 


M12*M22 

M21*M22 

M12*M32 

M31*M32 

M12*M42 

M41*M42 

M22*M32 

M3 1  *M3  2 

M22  *M4  2 

M4 1  *M4  2 

M32*M42 

M4 1 *M4  2 

*M22+M12*M21 ) 

*M22+M12*M21 ) 

*M42+M12*M41  ) 

*M42  +  M22*M41 ) 


•  • 


L23»DPPV*PPV+DPVV*PVV+DPPR*PPR+DPRR*PRR+DPPY*PPY+DPYY*PYY 
+  DVVR*VVR+DVRR*VRR+DVVY* VVY+DVYY* VYY+DRRY*RRY+DRYY*  RYY 
+  DPVR*  PVR+DPVY*  PVY  +  DPRY*  PRY+DVRY* VRY+DPPP*  PPP+DVVV* VVV 
+DRRR*RRR+DYYY* YYY 


L24-  DPPV*M12*M12*M22  +  DPVV*M12*M22*M22 
+  DPRR*M12*M32  *M32  +  DPPY*Ml2*Ml2*M42 
+  DWR*M22*M22*M32  +  DVRR*M22*M32*M32 
+  DVYY*M22*M42*M42  +  DRRY*M32*M32*M42 
+  DPVR*M12*M22*M32  +  DPVY*Ml2*M22*M42 
+  DVRY*M22*M42*M32  +  DPPP*Ml 2 *Ml 2*M1 2 
+  DRRR*M32*M32*M32  +  DYYY*M42*M42*M42 


DPPR*M12*M12*M32 
DPYY*M12*M42*M42 
DVVY*M22*M22*M42 
DRYY*M32*M42*M42 
DPRY*M12*M42*M32 
DVVV*M22*M22  *M22 


L31-L21 

L32-L22 

L33-L23 

L34-L24 


L21-L21*BB1*U*U 

L22-L22*BB1*U*U 

L23-L23*BB1*U*U 

L24«L24*BBl*U*U 

L31-L31*BB2*U*U 

L32-L32*BB2*U*U 


•  • 


non  non  ono 


L33-L33*BB2*U*U 

L34-L34*BB2*U*U 

L41  — 0.5*Mll*MllMM21+U*Ml 1/3.0) 

L42— M11*(M12*M21+0.5*M11*M22+0.5*U*M12*M12) 
L43— M12*(Mll*M22+0.5*M12*M21+0. 5*U*M11*M12 ) 
L44  — 0.5*M12*M12*(M22+U*M1 2/3.0) 

Invert  Transformation  Matrix 

DO  2  1-1,4 
DO  3  J-1,4 

TINV ( I , J ) -0 . 0  UuOt 

TSAVE ( I , J ) -T ( I ,  J  ) 

3  CONTINUE 
2  CONTINUE 

CALL  DLUD (4,4, TSAVE , 4 , TLUD , IVLUD ) 

DO  4  1-1,4 

IF  ( IVLUD ( I ) . EQ . 0 )  STOP  3003 

4  CONTINUE 

CALL  DILU (4,4, TLUD , IVLUD , SVLUD ) 

DO  8  1-1,4 
DO  9  J-1,4 

TINV ( I , J ) -TLUD( I , J ) 

9  CONTINUE 

8  CONTINUE 


Check  Inv ( T ) *A*T 


IMULT-2 

IF  (IMULT.EQ.l)  CALL  MULT( TINV, ASAVE , T  ) 
IF  (IMULT.EQ.0)  STOP 


C 


C 


Definition  of  Nij 

Nl 1— TINV ( 1,1) 
N21-TINV (2,1) 
N31-TINV (3,1) 

N41— TINV( 4,1 ) 
N12-TINV( 1,2) 
N22-TINV( 2,2) 
N32-TINV ( 3,2) 
N42-TINV ( 4,2) 

Nl 3-TINV ( 1,3) 
N23-TINV( 2,3) 
N33-TINV( 3,3) 

N4  3-TINV (4,3) 
N14-TINV( 1,4) 
N24«TINV( 2,4) 
N34-TINV ( 3,4) 
N44-TINV ( 4,4) 


Rll-N12*L21+N13*L31+N14*L41 
R12-N12  *L22+Nl 3*L32+N1 4  *L4  2 
R13-N12*L23+N13*L33+N14*L43 
R14-N12*L24+N13*L34+N14*L44 
R21-N22*L21+N23*L31+N24*L41 
R22-N22*L22+N23*L32+N24*L42 
R23-N22*L23+N23*L33+N24*L43 
R24-N22*L24+,J2  3*L34+N2  4*L4  4 


•  • 
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Evaluate  Alpha'  and  Oaega' 

DINOO.OOl 
XDR  -XD+DINC 
XDL  -XD-DINC 
XD  -XDR 
PSI  -0 . DO 

Y  -0 . DO 

V  -0.D0 
DPSI-KlP 
DV  -K2P 
DR  -K3 

DY  -R1*XD/(XD**2+Y**2) 

A( 1 , 1 )-0 . 0D0 
A( 1 , 2 )-0 . 0D0 
A( 1 , 3 ) -1 . 0D0 
A( 1 , 4 ) -0 . 0D0 

A( 2 , 1 )-  BB1 *U*U*DPSI 

A( 2 , 2 ) — AAl 1 *U+BB1 *U*U*DV 
A( 2 , 3 )-AAl2*U+BBl*U*U*DR 
A( 2 , 4 )-  BBl *U*U*DY 

A( 3 , 1 )-  BB2*U*U*DPSI 

A(  3 , 2 )«AA21*U+BB2*U*U*DV 
A( 3 , 3 )-AA22*U+BB2*U*U*DR 
A(  3 , 4  )  -  BB2*U*U*DY 

A(  4 , 1 )-U*DCOS( PSI )-V*DSIN( PSI ) 

A( 4 , 2 ) -  DCOS(PSI) 

A( 4 , 3 )-0 . 0D0 
A( 4 , 4 ) -0 . ODO 

CALL  RG( 4 , 4 , A,WR,WI , 0 , ZZZ , IVl.FVl, I  ERR) 
CALL  DSTABL( DEOS , WR , WI , FREQ ) 

ALPHR-DEOS 

OMEGR-FREQ 

XD-XDL 
PSI  -0 . DO 

Y  -0 . DO 

V  — 0 . DO 
DPSI-KlP 
DV  -K2P 
DR  -K3 

DY  -K1*XD/(XD**2+Y**2) 

A(  1 , 1 )  — 0 . ODO 
A(  1 , 2  )  -0 . ODO 
A( 1 , 3 )-l . ODO 
A( 1 , 4 )— 0 . ODO 

A(  2 , 1 )  —  BBl*U*U*DPSI 

A( 2 , 2 )-AAll*U+BBl*U*U*DV 
A( 2, 3 )-AAl2*U+BBl*U*U*DR 
A( 2 , 4 )-  BBl *U*U*DY 

A( 3 , 1 )-  BB2*U*U*DPSI 

A( 3 , 2 )— AA21*U+BB2*U*U*DV 
A( 3, 3 )-AA22*U+BB2*U*U*DR 
A( 3 , 4 )-  BB2*U*U*DY 

A( 4, 1 )-U*DCOS{ PSI )-V*DSIN( PSI ) 

A( 4 , 2 )»  DCOS(PSI) 

A( 4 , 3 ) -0 . ODO 


on  non 


cm  wm  i ,  <f,  A,  ww,wi ,  u,  ill  ,  1  n  ,tvi,  man 

CALL  DSTABL ( DEOS , WR ,  WI , FREQ ) 

ALPHL-DEOS 

OMEGL-FREQ 

DALPHA-( ALPHR-ALPHL )/( XDR-XDL ) 

DOMEGA- ( OMEGR-OMEGL ) / ( XDR-XDL ) 

Evaluation  of  Hopf  Bifurcation  Coefficients 


COEFl-3 .0*R11+R13+R22+3.0*R24 
COEF2-3  .0*R21+R23-R12-3.0*R14 
AMU 2  — COEF1/(8.0*DALPHA) 

BETA2-0 . 25*COEFl 

TAU2  — ( C0EF2-D0MEGA*C0EFl/DALPHA ) / ( 8 . 0  *OMEGAO ) 

PERD  -2.0*3. 141 59 2 7/OMEGA 0 

PER  — PERD*U/L 

X-XD/L 

K4-K1/XD 

WRITE  (11,*)  TC , COEFl 
WRITE  (12,*)  TC , PER 
WRITE  (13,*)  TC , DALPHA 
TC-TC+0 .015 
40  CONTINUE 

TL-TL+0 . 02 
50  CONTINUE 


1001  FORMAT  ( '  ENTER  TIME  CONSTANT' ) 

1005  FORMAT  ('  ENTER  TIME  LAG') 

1006  FORMAT  ('  ENTER  DO') 

2002  FORMAT  (5E15.5) 

END 

C— ........................................... 

SUBROUTINE  DSOMEG ( I JK , WR , WI , OMEGA , CHECK ) 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
DIMENSION  WR( 4 ) ,WI ( 4 ) 

CHECK— 1.0D+25 
DO  1  1-1,4 

IF  ( WR (I).LT.CHECK)  GO  TO  1 
CHECK-WR ( I ) 

IJ-I 

1  CONTINUE 

OMEGA-DABS ( WI ( I J ) ) 

IF  (WI(IJ) .GT.0.D0)  I JK— I J 
IF  (WI(IJ) .LT.0.D0)  IJK-IJ-1 
RETURN 
END 


SUBROUTINE  DSTABL ( DEOS , WR , WI , OMEGA) 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
DIMENSION  WR( 4 ) , WI ( 4 ) 

DEOS— 1 .0D+20 
DO  1  1-1,4 

IF  (WR( I ) .LT.DEOS)  GO  TO  1 
DEOS— WR ( I ) 

IJ-I 

1  CONTINUE 
OMEGA-WI ( IJ ) 

OMEGA-DABS ( OMEGA ) 
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RETURN 

END 


SUBROUTINE  NORMAL ( T ) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DIMENSION  T(4,4) ,TNOR(4,4) 

CFAC-T( 1 , 1 ) **2+T( 1 , 2 ) **2 
IT  (DABS(CrAC) .LE. (l.D-10) )  STOP  4001 
TNOR( 1 , 1 )-l .DO 

TNOR( 2 , 1 )-(T( 1 , 1 ) *T( 2 , 1 ) +T( 2 , 2 ) *T( 1 , 2 ) )/CFAC 
TN0R(3,1)-(T(1,1)*T(3,1)+T(3,2)*T{1,2) )/CFAC 
TNOR( 4,1)»(T(1,1)*T(4,1)+T(4,2)*T(1,2) )/CFAC 
TNOR( 1 , 2 ) -0 . DO 

TNOR(2,2)-(T(2,2)*T(l,l)-T(2,l)*T(l,2) )/CFAC 
TNOR(3,2)-(T(3,2)*T(l,l)-T(3,l)*T(l,2) )/CFAC 
TNOR(4,2)-(T(4,2)*T(l,l)-T(4,l)*T(l,2) )/CFAC 
DO  1  1-1,4 
DO  2  J-1,2 

T ( I , J )-TNOR ( I , J ) 

2  CONTINUE 
1  CONTINUE 
RETURN 
END 


SUBROUTINE  MULT ( T I NV , A , T ) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DIMENSION  TINV ( 4,4),A(4,4),T(4,4),Al(4,4),A2(4,4) 
DO  1  1-1,4 
DO  2  J-1,4 
Al ( I , J )— 0 .DO 
A2  ( I , J ) -0 . DO 

2  CONTINUE 
1  CONTINUE 

DO  3  1-1,4 
DO  4  J-1,4 
DO  5  K-1,4 

Al  (  I , J ) — A( I ,K)*T(K,J)  +A1  ( I  ,  J  ) 

5  CONTINUE 
4  CONTINUE 

3  CONTINUE 
DO  6  1-1,4 

DO  7  J-1,4 
DO  8  K-1,4 

A2 ( I , J ) -TINV ( I , K ) *Al ( K , J ) +A2 ( I , J ) 

8  CONTINUE 

7  CONTINUE 

6  CONTINUE 

DO  11  1-1,4 

WRITE  (*,101)  (A( I , J) , J-1,4) 

11  CONTINUE 

DO  12  1-1,4 

WRITE  (*,101)  (T(I,J) , J-1,4) 

12  CONTINUE 
DO  10  1-1,4 

WRITE  (*,101)  (A2( I, J) , J-1,4) 

10  CONTINUE 

WRITE  (*,101)  A2 (1,1) 

RETURN 

101  FORMAT  (4E15.5) 

END 
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APPENDIX  B 


C  PROGRAM  SIMU. FOR  (VEHICLE  PATH  SIMULATION) 

REAL  K1 , K2 , K3 , L , NR , NV , NDRS , NDRB , I Z , MASS , 

&  NRDOT,NVDOT,!UP,X2P 

c 

OPEN  ( 20, FILE-' SIMl.DAT' , STATUS- ' NEW' ) 

C  OPEN  ( 30,riLE»'SIM2.DAT' , STATUS- ' NEW ' ) 

C 

WEIGHT-435.0 
12  -45.0 

L  -7.3 
RHO  -1.94 
G  -32.2 
XG  -0.0104 
MASS  -WEIGHT/G 
U  -5.0 
RATIO  -0.0 
DO  -0.4 
XG  -XG/L 

MASS  -MASS/(0.5*RHO*L**3) 

12  -IZ/(0.5*RHO*L**5) 

YRDOT  —0.00178 
YVDOT  —0.03430 
YR  -+0.01187 
YV  —0.03896 
YDRS  -+0.02345 
YDRB  -+0.02345 
NRDOT  —0.00047 
NVDOT  —0.00178 
NR  —0.01022 
NV  —0.00769 
NDRS  —0.337*0.02345 
NDRB  -+0.283*0.02345 
C 

DH  — ( IZ-NRDOT ) * ( MASS-YVDOT ) - 
&  ( MASS*XG-YRDOT ) * ( MASS*XG-NVDOT ) 

AAl 1- ( ( IZ-NRDOT) *YV- ( MASS*XG-YRDOT ) *NV)/DH 
AA12- (  (  IZ-NRDOT)  M-MASS+YR)- 
&  ( MASS*XG-YRDOT ) * ( -MASS*XG+NR ) )/DH 

AA21- ( ( MASS-YVDOT ) *NV- ( MASS*XG-NVDOT ) * YV ) /DH 
AA22- ( ( MASS-YVDOT) *( -MASS *XG+NR)- 
&  ( MASS*XG-NVDOT ) * ( -MAS5+YR ) )/DH 

BBl 1- ( ( IZ-NRDOT) *YDRS- ( MASS*XG- YRDOT ) *NDRS )/DH 
BB12- ( ( IZ-NRDOT ) *YDRB- ( MASS *XG- YRDOT ) *NDRB ) /DH 
BB21- ( ( MASS-YVDOT ) *NDRS- ( MASS*XG-NVDOT ) * YDRS ) /DH 
BB22- ( ( MASS-YVDOT ) *NDRB- ( MASS*XG-NVDOT ) * YDRB ) /DH 
C 

BBl-BBl 1+RATI0*BB12 
BB2-BB21+RATIO*BB22 
C 

WRITE  (*,*)  'ENTER  TC , TL ,  D' 
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READ  (*,*)  TC,TL,D 

POLE-l.O/TC 

ADl-3.0*POLE 

AD2  - 3.0*  POLE*  *  2 

AD3-POLE**3 

Al-BBl 

B1-BB2 

C1--AD1- ( AA11+AA22 ) 
A2-(BBl*AA22-BB2*AA12) 

B2— ( BB2*AAll-BBl *AA21 ) 

Kl-AD3/( BB2*AAll-BBl*AA21 ) 

C2-AD2- ( AA1 1 *AA22-AAl 2 *AA2 1 ) +  BB2 *Kl 
K2-(C1*B2-C2*B1 )/( Al*B2-A2*Bl ) 

K3— ( C2*Al-Cl*A2 )/( Al*B2-A2*Bl ) 

TIME-100.0 
DT-0.1 
N-TIME/DT 
I Y-TL/DT 

PSI-  0.0 
V-  0.0 
R-  0.0 
Y  -  0.5 
T  -  0.0 
I COUNT- 0 
YCON-Y 

DO  10,  I— 1 , N 
T-  T+DT 

PSIDOT-R 

VDOT  — AAl 1 *V  +  AA12*R  +  BBl *DEL 
RDOT  -AA21 *V  +  AA22*R  +  BB2*DEL 
XDOT  —COS ( PSI ) -V*SIN ( PSI ) 

YDOT  *SIN( PSI ) +V*COS ( PSI ) 

PSI-PSI+ ( PSIDOT*DT ) 

V— V+ ( VDOT*  DT ) 

R— R+ ( RDOT*DT ) 

X— X+ ( XDOT*DT ) 

Y-Y+ ( YDOT*DT ) 

I COUNT— I COUNT+ 1 

IF  ( ICOUNT.NE.IY)  GO  TO  11 

YCON-Y 

I COUNT- 0 

PSIC— ATAN  (  YCON/D  ) 

DEL-Kl* (PSI-PSIC )+K2*V+K3*R 
DEL-DO  *TANH ( DEL/DO ) 

WRITE( 20 , * )  T , Y 
WRITE( 30 , * )  T, Y 
CONTINUE 
END 
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APPENDIX  C 


PROGRAM  THESNEWT . FOR  (TWO  POINT  TIME  DELAY  ANALYSIS) 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DOUBLE  PRECISION  Kl , K2 , K3 , L , NR , NV, NDRS , NDRB , IZ , MASS , 
&  NRDOT , NVDOT 

DIMENSION  A( 4 , 4  ) , FVl (4) , IV1 ( 4 ) ,  ZZZ( 4 , 4  ) ,WR( 4 ) , WI ( 4 ) 
DIMENSION  DIST( 10,10) 

OPEN  ( 11 , FILE- 'TCN8.DAT' , STATUS- ' NEW' ) 


Vehicle  Parameters: 
U-  5.0 
RATIO-  0.0 
WEIGHT-435.0 


IZ 

L 

RHO 

G 

XG 

MASS 

YRDOT 

YVDOT 

YR 

YV 

YDRS 

YDRB 

NRDOT 

NVDOT 

NR 

NV 

NDRS 

NDRB 


-45.0 

-7.3 

-1.94 

-32.2 

-0.0104 

-WEIGHT/G 

—  0.00178* 

—  0.03430* 
-+0. 01187* 

—  0.03896* 
-+0. 02345* 
-+0. 02345* 

—  0.00047* 

—  0.00178* 

—  0.01022* 
—  0.00769* 
—0.337*0. 
-+0.283*0. 


*0 . 5*RHO*L**4 

*0. 5*RHO*L**3 

*0 . 5*RHO*L**3 

*0 . 5*RHO*L**2 

*0 . 5*RHO*L**2 

*0 . 5*RHO*L**2 

*0. 5*RHO*L**5 

*0 . 5*RHO*L**4 

*0 . 5*RHO*L**4 

*0. 5*RHO*L**3 

. 02345*0. 5*RHO*L**3 

.02345*0. 5*RHO*L**3 


DH  -(IZ-NRDOT)* (MASS-YVDOT )- 
t  ( MASS*XG- YRDOT ) * ( MASS*XG-NVDOT ) 

AAll-( ( I Z-NRDOT ) *YV- ( MASS*XG-YRDOT ) *NV ) /DH 
AA12- ( ( IZ-NRDOT ) * ( -MASS+YR ) - 
i  ( MASS*XG-YRDOT ) * ( -MASS*XG+NR ) )/DH 

AA21- ( ( MASS-YVDOT ) *NV- ( MASS*XG-NVDOT ) * YV ) /DH 
AA22-( ( MASS-YVDOT)* (-MASS *XG+NR)- 
i  ( MASS  *XG-NVDOT ) * ( -MASS+YR ) ) /DH 

BBll- ( ( IZ-NRDOT ) * YDRS- ( MAS S*XG- YRDOT ) *NDRS ) /DH 
BB12- ( ( IZ-NRDOT ) * YDRB- ( MAS S*XG- YRDOT ) *NDRB ) /DH 
BB21- ( ( MASS-YVDOT ) *NDRS- ( MASS*XG-NVDOT ) * YDRS ) /DH 
BB22- ( ( MASS-YVDOT ) *NDRB- ( MASS*XG-NVDOT ) * YDRB ) /DH 

BBl— BBl 1+RATI0*BB12 
BB2-BB21+RATIO*BB22 


WRITE  (*,1005) 
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READ  ( * , * ) 


TL 


WRITE  (*,1006) 
READ  (  *  ,  *  )  TC 

EPS-1. D-10 
I THAX- 100000 


TL— 0 . 01 
DO  4  1-1,100 
TLD  -  TL  *  L/U 

DO  1  J-1,100 
TCD-TC*L/U 
POLEC— 1 . 0/TCD 
ADl-3 . 0*POLEC 
AD2-3 . 0*POLEC**2 
AD3-POLEC**3 
Al— BBl *U*U 
B1-BB2*U*U 

Cl— AD1-(AA11+AA22  )*U 
A2-(BB1*AA22-BB2*AA12 )*U**3 
B2«(BB2*AA11-BB1*AA21 )*U**3 
Kl-AD3/( (BB2*AA11-BB1*AA21)*U**3) 

C2— AD2-(AA11*AA22-AA12*AA21 ) *U**2+BB2*U*U*Kl 
K2-(C1*B2-C2*B1 )/( Al*B2-A2*Bl ) 
K3«(C2*A1-C1*A2)/(A1*B2-A2*B1 > 


& 


Al-  ( AAl 1 *BB2-AA2 l*BBl)*U*U*U*Kl 

A2-  (AA11*AA22-AA12*AA21 ) *U*U+ ( AAl 1 *BB2-AA21 *BBl ) *U*U*U*K3+ 
( AA22*BB1-AA12*BB2 ) *U*U*U*K2-BB2*U*U*K1 
A3— ( AA11+AA22 )*U-(BBl*K2+BB2*K3 ) *U*U 
BO-  ( AAl 1 *BB2-AA2 1 *BB1 )*U*U*U*U*Kl 
Bl— (AA12*BB2-AA22*BB1 ) *U*U*U*Kl-BB2*U*U*U*Kl 
B2— BBl*U*U*Kl 


•  • 


AQ-B2*A3-B1 

BQ-B1*A2-B2*A1-B0*A3 

CQ-B0*A1 

DET-BQ**2-4 . D0*AQ*CQ 

IF  (DET.LT.0.D0)  GO  TO  4 

ZTl- ( -BQ+DSQRT( DET ) )/(2.0D0*AQ) 

ZT2- ( -BQ-DSQRT ( DET ) )/(2.0D0*AQ) 

IF  (ZT1.GE.0.D0)  OM«DSQRT(ZTl ) 

IF  (ZT2.GE.0.D0)  OM-DSQRT( ZT2 ) 

ALPHAl-AQ 

ALPHA2-BQ 

ALPHA3-CQ 

BETA1— B2 

BETA2-B0+B2*A2-B1*A3 

BETA3-B1*A1-B0*A2 


OMOLD-OM  ^ 

3  BETA-BETAl *OMQLD*  *  5+BETA2  *OMOLD*  *  3  +  BETA3  *OMOLD 
ALPHA-ALPHAl*OMOLD**4+ALPHA2*OMOLD**2+ALPHA3 
BEP-5 . DO*BETAl*OMOLD**4  +  3 . DO*BETA2*OMOLD*  *  2  + BETA 3 
ALP-4 .DO*ALPHAl*OMOLD**3+2 . D0*ALPHA2 *OMOLD 
F 1-BETA* DS IN ( 2 . DO*OMOLD*TLD ) -2 . D0*BETA*DSIN( OMOLD*TLD ) 

F2 -ALPHA *DCOS ( 2 . DO*OMOLD*TLD ) -2 . D0*ALPHA*DCOS ( OMOLD*TLD )  • 
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o  o 


r-  n+p2 

FPl«BEP*DSIN( 2 .DO*OMOLD*TLD)+BETA*2 . DO*TLD*DCOS ( 2 . DO*OMOLD*TLD ) 

PP2— 2 . DO*BEP*DSIN( OMOLD*TLD ) -2 . DO*BETA*TLD*DCOS ( OMOLD*TLD ) 

FP3-ALP*DC0S ( 2 . DO*OMOLD*TLD ) -ALPHA*2 . DO*TLD*DSIN( 2 . DO*OMOLD*TLD ) 

PP4— 2 . DO*ALP*DCOS ( OMOLD*TLD ) +2 . DO  *TLD*ALPHA*DSIN ( OMOLD*  TLD ) 

FPR1ME-  FP1+FP2+FP3+FP4 

ir  (FPFIME.EQ.O.DO)  STOP  1112 

IF  (F.EQ.O.DO)  GO  TO  2 

OMNEW-OMOLD-F/FPRIME 

OMDI F-DABS ( OMNEW-OMOLD ) 

IF  (OMDIF.LE.EPS)  GO  TO  2 

OMOLD-OMNEW 

IT-IT+1 

IF  ( IT.GT. ITNAX)  STOP  1111 
GO  TO  3 

2  OM-OMNEW 

XDNUMl-DSQRT< ( BO-B2*OM* *2 ) * *2+ ( Bl *OM ) * *2  ) 

XDNUM2-DSQRT ( ( 2 . DO*DCOS ( OM*TLD ) -DCOS ( 2 . DO  *OM*TLD ) ) *  *2  + 

&  ( DSIN ( 2 . DO*OM*TLD ) -2 . DO  *DS IN ( OM*TLD ) )**2) 

XDNUM-XDNUMl *XDNUM2 

XDDEN-OM*DSQRT( ( A1-A3 *OM* * 2 ) *  * 2+ ( A2 *OM-OM* *  3 ) * * 2 ) 

XD-XDNUM/XDDEN 

X-XD/L 

C  DIST(I,J)-X 

WRITEdl,*)  X ,  TL 
TC-TC  +  0.015 
1  CONTINUE 
TL-TL  +  0.02 
4  CONTINUE 

C  WRITE(11, 20)  ( (DIST(I,J) ,J«1, 10) ,1-1,10) 

20  FORMAT( 10 ( 2X , F4 . 2  )  ) 

1005  FORMAT  ('  ENTER  TIME  LAG') 

1006  FORMAT  ('  ENTER  TIME  CONSTANT') 

END 
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